clc,clear;
close all;
%导入相关数据
data=xlsread("random_cities_coordinates.xlsx");
% 城市坐标
cities=data;
% 城市数量
nCities = length(cities);
%计算距离矩阵
dist = zeros(nCities, nCities);
for i = 1:nCities
    for j = 1:nCities
        dist(i, j) = sqrt(sum((cities(i, :) - cities(j, :)).^2));
    end
end
maxiter=1000; %最大迭代次数
startcity=5;    %初始城市
[bestx,bestf,bd]=sa(nCities,dist,maxiter,startcity);
disp("Best Solvution");
disp(bestx);
disp("Best Function");
disp(bestf);
tusi(bd,nCities,cities,bestx);

%%  sa
function [best_r, best_d,bd] = sa(ncities, dist,maxiter,startcity)
    %参数设置
    intemp=1000;    %初始温度
    rate=0.99;  %冷却率
    
    % 初始化随机路径
    recities = setdiff(1:ncities,startcity);
    croute = [startcity, recities(randperm(length(recities)))];
    cdistance = caldistance(croute, dist);
    
    % 设置初始最佳路径
    best_r = croute;
    best_d = cdistance;
    
    for i = 1:maxiter
        % 量子退火过程
        temp = intemp;
        while temp>0.01
            % 随机扰动生成新路径
            new_route = perroute(croute);
            new_distance = caldistance(new_route, dist);
            
            % 根据能量差决定是否接受新路径
            delta_energy = new_distance - cdistance;
            
            if delta_energy < 0 || exp(-delta_energy / temp) > rand
                croute = new_route;
                cdistance = new_distance;
            end
            
            % 更新最佳路径
            if cdistance < best_d
                best_r = croute;
                best_d = cdistance;
            end
            
            % 降低温度
            temp = temp * rate;
        end
        bd(i)=best_d;
    end
    %最优路径
    best_r=[best_r,startcity];
end

function total_distance = caldistance(route, distance_matrix)
    % 计算给定路径的总距离
    total_distance = 0;
    for i = 1:length(route)-1
        total_distance = total_distance + distance_matrix(route(i), route(i+1));
    end
    total_distance = total_distance + distance_matrix(route(end), route(1)); % 返回起点
end

function new_route = perroute(route)
    % 随机交换路径中的两个城市
    idx = randperm(length(route) - 1, 2) + 1; % 从第二个城市开始随机选择
    new_route = route;
    new_route(idx(1)) = route(idx(2));
    new_route(idx(2)) = route(idx(1));
end
%% 作图

function ktu=tusi(l,nCities,cities,bP)
l=[l(2:end)];
n=length(l);
count=1:n;
x=zeros(nCities+1,1);
y=zeros(nCities+1,1);
for i=1:nCities
    x(i)=[cities(bP(i),1)];
    y(i)=[cities(bP(i),2)];
end
x(i+1)=cities(bP(1),1);
y(i+1)=cities(bP(1),2);

figure(1)
plot(count,l);
xlabel('迭代次数');
ylabel('最优值');
title('模拟迭代图');
hold on
figure(2)
plot(x,y,'-ro', 'LineWidth', 2, 'MarkerEdgeColor', 'k', 'MarkerFaceColor', 'r', 'MarkerSize', 10);
xlabel('X 坐标');
ylabel('Y 坐标');
title('路径图');
grid on; % 打开网格
end
